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We present an overview of our project of large-scale simulations with dynamical overlap fermions. 
The first production run in two-flavor QCD is on-going using the Iwasaki gauge action on a 
16 3 x 32 lattice at the lattice spacing of 0.12 fm with six sea quark masses down to m s ,ph ys /6, 
where m S]P hy S is the physical strange quark mass. We briefly introduce our choice of the lattice 
action and simulation algorithm, and describe the present status of the production run. Preliminary 
results on the light meson masses and the static quark potential are also reported. 
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1. Introduction 

The JLQCD Collaboration has pursued detailed studies of non-perturbative properties of QCD 
through numerical simulations on the supercomputer system at KEK. Our previous unquenched 
simulations M, ||] however were restricted to rather heavy sea quark masses > m i)P hy S /2, which 
make the chiral extrapolation difficult to control. The use of Wilson-type quarks also limits the 
applicability of our simulation to quantities, such as Bk, for which the explicit chiral symmetry 
breaking causes complicated operator mixing. 

This year, the supercomputer system at KEK was replaced by a multiple system of Hitachi 
SRI 1000 and IBM Blue Gene/L with a total peak speed of about 60 TFLOPS. On this new system, 
we carry out large-scale simulations including dynamical overlap fermions with sea quark masses 
down to m.y phys/4 or smaller on reasonably fine {a < 0.125 fm) and large lattices (L>2 fm). In this 
article, we present an overview of our first production simulation in two-flavor QCD. 

2. Simulation method 

Our numerical simulations are carried out using the overlap quark action with the standard 
Wilson kernel Hw 

D(m se3 ) = (mo + + (m - y 5 sgn [H w (-m Q )] , (2. 1) 

where we set mo = 1.6. A major problem with this kernel is the appearance of (near-)zero modes 
of //w in simulations on relatively coarse lattices. This possibly spoils the locality of D and also 
makes simulations costly. From our preparatory study in quenched QCD we adopt the Iwasaki 
gauge action, with which the near-zero modes are suppressed compared to the standard plaquette 
glue and the localization range of D is sufficiently small. 

The use of the improved gauge action does not completely rule out the appearance of exact zero 
modes. In order to avoid the time consuming reflection/refraction procedure [Q], we introduce two- 
flavors of unphysical extra- Wilson fermions with a large negative mass —mo []5j |6|] and additional 
twisted mass ghosts ^ which produce the Boltzmann weight 

det [H w {-mof] 
det[# w (-m ) 2 + MT ' 

The factor in the numerator suppresses the zero modes during continuous updating of gauge fields, 
such as the hybrid Monte Carlo (HMC) algorithm. The denominator with appropriately chosen 
twisted mass /I cancels possibly large effects from high modes of Hw in the numerator. We note 
that these unphysical extra-fields have masses of cutoff order, and hence their effects vanish in the 
continuum limit. The suppression of zero modes in quenched and two-flavor QCD is demonstrated 
inRefs.[§,[7|]. 

While this procedure fixes the net topological charge Q during the HMC update, it does not 
forbid local topological fluctuations. It is expected that, in the infinite volume limit, properties 
of hadrons with their size of 0(1 fm) are insensitive to the global topological charge of gauge 
configurations, i.e. the effect due to the fixed topology is a finite size effect (FSE). We refer to 
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Ref.|7|] for more detailed discussions. We also note that our lattice action to fix the topology 
provides a convenient framework for studies in the £ -regime [gj. 

In our simulation program, we implement multiplications of D employing the low mode pre- 
conditioning: namely, we introduce a threshold An, in the spectrum of Hyj, and eigenvalues smaller 
than Ah are projected out in the evaluation of sgn[Z%]. For higher modes, we use the rational ap- 
proximation with the Zolotarev coefficients and the multi-shift CG algorithm for the inner-loop [Q]. 
For multiplications of D we employ a nested solver with the relaxed CG [fR]] for the outer-loop. 

In our implementation of HMC, we employ a combination of the Hasenbusch preconditioning 



[11] and the multiple time scale discretization of the molecular dynamics (MD) [12], which has 



been shown to be effective in simulations with Wilson-type fermions Q13|]. The determinant factor 
for overlap quarks is written as 



det [D(m sea ) 2 ] = det [D(m'f] det 



D(m s , 



\2l 



D(m'Y 



(2.3) 



where m' is the mass of the Hasenbusch preconditioner. Two pseudo-fermions, which we call PF1 
and PF2 in the following, are introduced for the first and second factors in the r.h.s. Then, PF2 

(PF2) 

is updated A^d times per unit trajectory length. The updates of PF1 and gauge fields are more 

(PF) (PF) (G) 

frequent by factors of /? MD and ^md ^md- Since we suppress zero modes by the extra- Wilson 
fermion, we switch off the reflection/refraction procedure in our MD evolution. 

In addition to the above algorithmic techniques, we also use an assembler code developed by 
IBM for multiplications of Dw on Blue Gene/L. This code makes the best use of the double FPU 
instructions, which process double precision complex number arithmetic effectively using two sets 
of floating-point registers. This assembler code is a factor of 3 faster than our naive Fortran code. 



Further details of our simulation algorithm are presented in Ref. [14]. 



3. Production run 



Our first production run with two-flavors of degenerate up and down quarks is being carried 
out with the above mentioned lattice action. The twisted mass /I = 0.2 is chosen from a preparatory 
study in quenched QCD [0]. We simulate j8 =2.30, which is expected to correspond to the lattice 
spacing about 0.125 fm, on a 16 3 x 32 lattice. Six sea quark masses listed in Table [l] are taken to 





m' 


/V MD 


MD 


MD 


HMC traj. 


Phmc 


time[min] 


#conf pot 


#COnfhad 


0.015 


0.2 


9 


4 


5 


2150 


0.89 


6.1 


151 


173 


0.025 


0.2 


8 


4 


5 


4320 


0.90 


4.7 


389 


410 


0.035 


0.4 


6 


5 


6 


4150 


0.74 


3.0 


411 


403 


0.050 


0.4 


6 


5 


6 


3500 


0.79 


2.6 


310 


310 


0.070 


0.4 


5 


5 


6 


3500 


0.81 


2.1 


307 


243 


0.100 


0.4 


5 


5 


6 


3590 


0.85 


2.0 


301 


319 



Table 1: Simulation parameters for configuration generation. Statistics for the measurements of static po- 
tential (#conf P ot) and hadron correlators (#confh a d) are also listed. 
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explore a range from m s ^ y& down to m SjP h ys /6. All results in this article are obtained in the trivial 
topological sector £2 = 0. 




200 400 
HMC traj 



200 400 
HMC traj 



Figure 1: Average (left panel) and maximum value 
(right panel) of MD force at each trajectory. Data at 
?«sea = 0.015 are plotted. 



In all simulations, we fix Ath = 0.045, with 
which several eigenvectors are projected out in 
the calculation of sgn[//w]- The number of poles 
in the Zolotarev approximation is set to 10 lead- 
ing to an accuracy of |sgn[//w] 2 — 1 1 — 10~( 7 ~ 8 ) 
throughout the simulations. The stopping condi- 
tion of the overlap solver is chosen so that the re- 
versibility in the gauge field and the Hamiltonian 
is satisfied to a level comparable to our previous 
simulations. 

We are generating gauge configurations on 
Blue Gene/L and store them on disks every 10 
trajectories. The HMC trajectory length is fixed 

to 0.5. Hadron correlators and the static quark potential are measured on both of SRI 1000 and 
Blue Gene/L using the stored configurations. Current statistics are summarized in Table ||. 

Figure [T] shows an example of the time history of the MD force from gauge fields, PF1 and PF2. 
With our choice of m', these forces are well separated from each other. This situation suitable for 
the application of the multiple time scale discretization is confirmed also at other sea quark masses. 
In the same figure, we also plot the force from the extra-fields in Eq. (2.2). Its maximum value 
becomes as large as that for the gauge field probably due to small eigenvalues of //\y- Therefore, 
we update the extra-fermions in the inner most MD step. 

As summarized in Table [j], we achieve the so, — , — , — , — , — , — , — , — , — , — , — , — , 
acceptance rate Phmc of 74 % or higher with 
our choice of the MD discretization parameters. 
The CPU time per 1 trajectory on the whole ma- 
chine of Blue Gene/L (10 racks, 57.3 TFLOPS) 
is also listed in the table. It takes about one 
month to accumulate 2000 trajectories at six sea 
quark masses. We have recently switched to a 
five dimensional solver proposed in Ref.JTBj], which 
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Figure 2: Autocorrelation time for plaquette and it- 
eration count of the overlap solver Ni nv . 



reduces the CPU time by roughly a factor of 2 
for our simulation parameters. 

In Fig. ||, we plot the sea quark mass depen- 
dence of the autocorrelation time for two quan- 
tities: Tpi a q for the plaquette, and T; nv for the iteration count of the overlap solver N[ nv . We find that 
Tpiaq is small and has a mild mass dependence, probably because it is a local quantity. On the other 
hand, N- mv is expected to be sensitive to the low modes of D. In fact, Ti nv increases rapidly towards 
the chiral limit. 

We find that T; nv « 60 at our lightest sea quark mass. This implies that statistics of 0(10,000) 
trajectories are needed for precise calculation of hadronic observables at such light quark masses. 
In the following, we present preliminary results with limited statistics listed in Table [j]. 
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Figure 3: Plot of (V(r)— Vft)ro as a function of r/ro (left figure), and chiral extrapolation of lattice spacing 
determined from ro (right figure). 



4. Static quark potential 



We calculate the static quark potential V{r) from the smeared Wilson loops. The Sommer 
scale ro fll6[ ] is extracted from a parametrization V (r) = Vq — cc/r + a r. 

In Fig. H| we plot V (r) rescaled by ro at all sea quark masses. While we simulate a wide range 
of the sea quark mass, all data form a universal curve. We note that our statistics are not sufficiently 
high to observe a clear dependence of a on the sea quark mass. The chiral extrapolation of the 
lattice spacing determined from ro is shown in the same figure. From a simple linear fit and an 
input ro = 0.49 fm, we obtain a = 0.1 199(14) fm in the chiral limit, which is close to our target 
value 0.125 fm. 



5. Meson masses and decay constant 

Meson masses are extracted from correlation functions measured with smeared source and 
local sink operators. Figure |] shows the chiral extrapolation for the pseudo-scalar and vector 
mesons made of degenerate valence quarks. With our statistics, we do not observe any significant 
deviation from simple linear fits 
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Figure 4: Chiral extrapolation of pseudo-scalar (left figure) and vector meson masses (right figure). 
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which give # 2 /dof < 1.0. Employing this linear fit, we obtain a=0. 1312(23) fm from the p meson 
mass input. This is consistent with the estimate from ro within 10% accuracy. 

Since we are exploring small sea quark mass regime with a fixed topology, we have to be 
careful about FSEs. These are however difficult to assess reliably without performing direct sim- 



ulations with different lattice volumes. Here we only note that an analytic formula for FSE [ |17[ ] 
combined with an analysis for the fixed topology Jl8| ] suggests that FSEs are not large 1 - 2%) 
even at our lightest sea quark mass. 

In order to check the size of the chiral sym- 
metry breaking, we calculate the following ratio 
motivated from the axial ward identity 

(V 4 A 4 P t ) 



Pawi 



(PPt) 



(5.2) 



where P and are the pseudo-scalar and ax- 
ial current operators. Fig. || shows the chiral ex- 
trapolation of Pawi with the following linear and 
quadratic forms: 

PAWI = «AWI + ^AWI »*sea (+ c AWI»4, a ) ■ ( 5 - 3 ) 




Figure 5: Chiral extrapolation of ratio Pawi defined 



in Eq. (5.2) 



The intercept ^awi gives a measure of the chiral symmetry breaking. We confirm that results from 
linear (0.0010(7)) and quadratic fits (-0.0007(9)) are consistent with zero. 

The slope b^wi provides a non-perturbative estimate of the renormalization constant Za = 
2/&awi- We obtain 1.394(20)(81), where the first error is statistical, and difference in Za between 
the linear and quadratic extrapolations is added as a systematic error. 

In Fig. ^, we plot our results for the pseudo- 
scalar decay constant together with our previ- 
ous estimates with the 0(a) -improved Wilson 
quarks [p. All results are renormalized non- 
perturbatively by Za = 2/&awi or Za in Ref.[ 19]. 
At heavy quark masses, we observe a reason- 
able consistency between our new and previous 
results with different discretizations. This sug- 
gests that these data are close to their continuum 
limit and have small scaling violations. 

At small sea quark masses, however, its quark Figure 6: Chiral behavior of / PS . Squares are previ- 
mass dependence is not smooth probably due to ous results with 0(a)-improved Wilson quarks [§. 
the limited statistics. As a result, chiral extrapo- 
lation with the chiral logarithmic term is very unstable. At this moment, it is difficult to make a 
definite conclusion on the consistency with chiral perturbation theory. 
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6. Summary 

In this article, we report on the JLQCD 's new project of large-scale simulations with dynamical 
overlap quarks. With careful choice of the lattice action and implementation of various algorithmic 
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techniques, we demonstrated that it is feasible to accumulate high statistics on reasonably fine 
(a ^ 0.125 fm) and large lattices (L«2 fm). 

Our simulations have just started. We are currently accumulating 0(10,000) trajectories for a 
precise determination of meson masses and decay constants, with which we are going to investigate 
the consistency with chiral perturbation theory. 

One of the most important subjects in the near future is the systematic effects due to the fixed 
topology. Simulations with Q / are underway to check Q dependence of physical observables. 
We are also planning to extend our simulations to three-flavor QCD. Pushing simulations toward 
larger volumes is also interesting future direction to explore the baryon sector as well as to check 
FSEs due to the small sea quark masses and the fixed topology. 

Numerical simulations are performed on Hitachi SRI 1000 and IBM System Blue Gene Solu- 
tion at High Energy Accelerator Research Organization (KEK) under a support of its Large Scale 
Simulation Program (No. 06-13). This work is supported in part by the Grant-in- Aid of the Ministry 
of Education (No. 13135204, 13135213, 13135216, 15540251, 16540228, 16740147, 16740156, 
17340066, 17540259, 17740171, 18034011, 18104005, 18340075, 18740167). 

References 

[1] S. Aoki et al. (JLQCD Collaboration), Phys. Rev. D68 (2003) 054502. 

[2] T. Ishikawa et al. (CP-PACS and JLQCD Collaborations), PoS (LAT2006) 181. 

[3] N. Yamada et al. (JLQCD Collaboration), PoS (LAT2006) 060. 

[4] Z. Fodor, S.D. Katz and K.K. Szabo, JHEP 0408 (2004) 003. 

[5] P.M. Vranas, arXiv:hep-lat/0001006; hep-lat/0606014. 

[6] H. Fukaya et al. (JLQCD Collaboration), arXiv:hep-lat/0607020. 

[7] S. Hashimoto et al. (JLQCD Collaboration), PoS (LAT2006) 052. 

[8] H. Fukaya et al. (JLQCD Collaboration), PoS (LAT2006) 050. 

[9] A. Frommer et al, Int. J. Mod. Phys. C6 (1995) 627. 
[10] N. Cundy et al., Comput. Phys. Commun. 165 (2005) 221. 
[11] M. Hasenbusch, Phys. Lett. B519 (2001) 177. 
[12] J.C. Sexton and D.H. Weingarten, Nucl. Phys. B380 (1992) 665. 

[13] M.J. Peardon and J. Sexton, Nucl. Phys. (Proc.Suppl.) 1 19 (2003) 985; A. Ali Khan et al. (QCDSF 
Collaboration), Phys. Lett. B564 (2003) 235; C. Urbach, K. Jansen, A. Shindler and U. Wenger, 
Comput. Phys. Commun. 174 (2006) 87. 

[14] H. Matsufuru et al, (JLQCD Collaboration), PoS (LAT2006) 031. 

[15] A. Borici, arXiv:hep-lat/99 12040; hep-lat/0402035; R.G. Edwards et al, PoS (LAT2005) 146. 

[16] R. Sommer, Nucl. Phys. B411 (1994) 839. 

[17] FC. Hansen and H. Leutwyler, Nucl. Phys. B350 (1991) 201. 

[18] R. Brower, S. Chandrasekharan, J.W. Negele and U.-J. Wiese, Phys. Lett. B560 (2003) 64. 
[19] M. Delia Morte et al. (ALPHA Collaboration), JHEP 0507 (2005) 007. 



7 



